Evidence for softening of first-order transition in 3D by quenched disorder 
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We study by extensive Monte Carlo simulations the effect of random bond dilution on the phase 
transition of the three-dimensional 4-state Potts model which is known to exhibit a strong first-order 
transition in the pure case. The phase diagram in the dilution-temperature plane is determined from 
the peaks of the susceptibility for sufficiently large system sizes. In the strongly disordered regime, 
numerical evidence for softening to a second-order transition induced by randomness is given. Here 
a large-scale finite-size scaling analysis, made difficult due to strong crossover effects presumably 
caused by the percolation fixed point, is performed. 



The influence of random, confining geometries on first- 

' order phase transitions has been the subject of exciting 
experimental studies in the past few years. The case of 
the isotropic to nematic transition olnCB liquid crystals 
confined into the pores of aerogels consisting of multiply 
connected internal cavities has been particularly exten- 
sively studied and led to spectacular results: The first- 

1 order transition of the corresponding bulk liquid crystal is 
drastically softened in the porous glass and becomes con- 
tinuous jlj, an effect which was not attributed to finite- 
size effects but rather to the influence of random disorder. 
The first attempt to reproduce such a softening sce- 

, nario using Monte Carlo (MC) simulations was reported 
by Uzelac et al. M who studied a three-dimensional (3D) 

, g-state Potts model with spin variables (taking q = 3 
and 4 states per spin) located inside the randomly con- 

, nected pores of an aerogel modeled by diffusion-limited 
cluster aggregation. Although in experimental studies [|| 

, commonly random fields or random uniaxial anisotropics 
are suggested to explain the softening of the transition, 

, the random disorder chosen in Ref. [|2| is coupled to the 

' energy density and thus more akin to bond-dilution. 
The qualitative effect of random bond disorder 
on second-order phase transitions is well understood 
through the Harris relevance criterion Q, and a beauti- 

' ful experimental confirmation was reported in a LEED 
investigation of a two-dimensional (2D) order-disorder 
transition Q . For systems with a first-order transition in 

, the pure case, randomness generically softens the transi- 

1 tion and, under certain circumstances, may even induce 
a second-order transition according to a picture first pro- 
posed by Imry and Wortis H . 

In 2D, the natural candidate for theoretical inves- 
tigations is the g-state Potts model, since this model 
is exactly known to exhibit regimes with first- and 
second-order transitions |Q, depending on the value of 
q. Many results were obtained in both regimes in the 
last ten years |^,^) , including approximate analytic treat- 
ments, MC simulations, transfer-matrix calculations, and 
high-temperature series expansions. Among others also 



quite intricate problems such as self-averaging and multi- 
fractality have recently been studied in some detail ||. 

In 3D, to date only the Ising model with site dilution 
has been studied extensively ju|. In accordance with 
the Harris criterion the presence of random disorder was 
found to modify the critical exponents to values close to 
v = 0.684(5), 7/1/ = 1.96(3), and (3/v = 0.519(8). Con- 
cerning the influence of random disorder on first-order 
transitions, even less is known in 3D. Apart from the 
exploratory work jllj of Uzelac et al. Q), only the site 
diluted 3-state Potts model, which in the pure 3D case 
has a very weak first-order transition [ ^2[ , has recently 
been studied via large-scale MC simulations E|. This 
study led to the conclusion that the critical exponent v 
governing the scaling behavior of the correlation length is 
compatible with that of the 3D site diluted Ising model, 
whereas the 77 exponent is definitely different. 

The purpose of this Letter is to present numerical ev- 
idence for softening of the transition when it is strongly 
of first order in the pure system, in order to be sensitive 
to disorder effects. The paradigm in 3D is the 4-state 
Potts model, since the correlation length at the tran- 
sition temperature of the unperturbed system is small 
enough (£ ~ 3 in lattice spacing units U) to allow 
simulations of significantly large systems. For the pure 
3D 5-state Potts model the first-order transition is al- 
ready too strong p3|. In the following we, therefore, 
consider the 4-state bond diluted Potts model on simple- 
cubic lattices of size V = L 3 with periodic boundary 
conditions. The Hamiltonian of the system with in- 
dependent, quenched random interactions is written as 
— (37i = j) Kij&<Ti.<jj, where the spins take the val- 
ues <7i = 1,...,4 and the sum goes over all nearest- 
neighbor pairs (i,j). The coupling strengths are allowed 
to take two different values Kij = K = J/ksT and 
with probabilities p and 1 — p, respectively. The or- 
der parameter for a given realization of the is de- 
fined by the majority orientation of the spins, m = (/i), 
where fi = {qp max -l)/(q — l) and p max is the maxi- 
mum value of the density of spins p a in the q = 4 pos- 



sible spin states. The thermal average over the MC it- 
erations is indicated by brackets (...), and the physical 
quantities are then averaged over disorder realizations, 
e.g., m = {(i). For each disorder realization the suscep- 
tibility is defined as usual via the fluctuation-dissipation 
theorem, X = KV({^ 2 ) - (n) 2 )- 

The present MC simulations consist of two parts. 
First, a scan of the dilution-temperature plane in order 
to determine the phase diagram, and second, a large- 
scale finite-size scaling (FSS) study at p — 0.56 up to 
L = 96 in the dilution regime exhibiting second-order 
transitions. The spin updates were performed with the 
cluster-flipping method E3| in the Swendsen-Wang for- 
mulation which turned out to be better behaved than the 
Wolff single-cluster version for high dilutions (small p), 
where small clusters connected by non-vanishing bonds 
are more likely to appear. We made sure that, by adjust- 
ing the length of the runs of cluster algorithms, at least 
10 tunnelling events between the two coexisting phases 
of the pure system (p — 1) were observed up to lattice 
size L = 16. To improve the accuracy, for weak dilu- 
tions (p between 1 and 0.68) where we obtained evidence 
for first-order transitions, we also used multicanonical al- 
gorithms For the determination of the maxima of 
observables we applied standard histogram reweighting 
techniques in order to extrapolate the results over a tem- 
perature range around the simulation point. 

At each probability p and for each realization of the 
random couplings, between 15 x 10 3 to 30 x 10 3 MC 
sweeps per spin were performed, resulting in at least 
250 (almost) independent measurements of the physical 
quantities for the largest lattice size considered. This 
turned out to be sufficient for reliable thermal averages. 
For the average over disorder realizations, between 2 000 
and 5 000 samples were generated. 

The phase diagram is determined from the locations 
of the maxima of the average susceptibility, i^ max , ob- 
tained for systems of increasing sizes up to L = 16. 
Theoretically we expect that the transition remains of 
first order in the regime of low impurity concentrations. 
For increasing concentrations a regime of second-order 
transitions should appear from a (tri-) critical concen- 
tration until the percolation threshold is reached where 
the transition vanishes altogether. After this percolation 
threshold, no ordered phase can exist at any finite tem- 
perature. Two points are known in the p — ksT / J plane: 
The transition temperature T t of the pure system fl4|| 
k B T t {p = I)/ J ~ 1.59076, and k B T t {p = p c )/J = 
at the percolation threshold p c w 0.2488. As can be 
inspected in Fig. [I], the temperatures of the suscepti- 
bility maxima for different lattice sizes are very stable 
and already for 10 3 spins an accurate transition line is 
obtained. Also shown is the result of a simple mean- 
field argument taking into account the average number 
of neighbors, kBT t (p)/zJ = const, where the constant 
is chosen such that T t of the pure system is reproduced. 



This leads to a simple linear approximation of the transi- 
tion line, kBT t (p)/ J = 1.59076 x p, which is surprisingly 
accurate over a significant range of p- values. 
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FIG. 1. Phase diagram of the bond diluted 4-state Potts 
model in 3D. The solid line is a spline interpolation to guide 
the eye, and the dashed line shows the transition line within 
a simple mean-field argument, T t (p) =pT t (l). 

By monitoring the FSS behavior of various thermo- 
dynamic quantities as well as the (pseudo-) dynamics of 
the update algorithm, we estimate the tricritical point 
to be located between p — 0.68 and 0.84: The shape of 
the energy probability densities as well as the Binder cu- 
mulants suggest a first-order transition at p — 0.84 and 
above while a clear second-order signal is observed at 
p = 0.68 and below. For the investigation of the criti- 
cal properties in the second-order regime we have chosen 
a dilution p = 0.56 where the corrections to scaling are 
seemingly relatively small, since the effective transition 
temperatures corresponding to the susceptibility maxima 
remain almost constant in the range of sizes 2 < L < 16 
used for the determination of the phase diagram. 

In order to convince ourselves that for p = 0.56 the 
transition is indeed of second order, let us first consider 
the average probability densities of the energy, P(e). In 
Fig. H their shapes close to K maK are depicted for vari- 
ous lattice sizes up to L = 96. We see that the system 
exhibits for small sizes two distinct peaks which clearly 
collapse into a single peak when one approaches the ther- 
modynamic limit. This is precisely what is expected at 
a second-order phase transition, while in the case of a 
first-order transition the double-peak structure should 
persist for all sizes and, in fact, should become even more 
pronounced when the system size increases. Physically, 
a two-peak structure would reveal the presence of two 
phases at the transition temperature, and this coexis- 
tence is the characteristic feature of a first-order tran- 
sition. More quantitatively we have confirmed that the 
effective interface tension a a derived from the usual re- 



2 



lation P m in/-Pmax ~ exp(— 2o- d vanishes in the 
infinite- volume limit; see the inset of Fig. |[ 

Our final goal is a quantitative characterization of the 
critical behavior by providing estimates for the critical 
exponents of the transition. To this end we have per- 
formed a standard FSS analysis at p = 0.56. As can be 
inspected in Fig. [| the corrections to asymptotic FSS 
for x ma x seem to become quite small above L = 30. 
The data are thus linearly fitted to a x L' y / 1 ' for L in the 
ranges L nlln > 30 to £ max = 96, and the resulting ex- 
ponents are collected in Table ||. Selecting the fits with 
the smallest chi-squared per degree of freedom, x 2 /d.o.f, 
we take as the final result the lines in bold face, e.g., 
7 /z/= 1.50 ±0.02. 
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FIG. 2. The probability density of the energy at p = 0.56 
for sizes L — 16, 20, 25, 30, 50, 64, and 96, reweighted to the 
temperature where the two peaks are of equal height. The 
inset shows the associated effective interface tension. 
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FIG. 3. Finite-size scaling benavior of the susceptibility, 
the magnetization and of dlnm/dK at A" max (the quantities 
have been shifted in the vertical direction for the sake of clar- 
ity). The scaling behavior for small lattice sizes is presumably 
governed by the percolation fixed point, and above a crossover 
length scale it reaches a new (random) fixed point. 



TABLE I. 


Linear fits for 


Xmax, (Ok In m)K m 


lx , and m_fc max • 


^min 


^max 


7/1/ 


x 2 /d.o.f 


35 


96 


1.500(14) 


0.044 


40 


96 


1.502(17) 


0.054 


50 


96 


1.506(27) 


0.065 


-^min 


^max 


1/u 


xVd.o.f 


35 


96 


1.362(13) 


1.011 


40 


96 


1.353(16) 


0.887 


50 


96 


1.330(25) 


0.419 






p/v 


xVd-o.f 


35 


96 


0.592(13) 


2.778 


40 


96 


0.608(15) 


2.145 


50 


96 


0.645(24) 


0.311 



The quantity &k In tu gives an estimation of the expo- 
nent v, (8k lnm)/f max oc L 1/>u . Here our analysis leads 
to an estimate ofl/v= 1.33 ±0.03 or v = 0.752 ±0.014, 
in agreement with the stability condition of the random 
fixed point (y > 2/D = 0.666...) and significantly dif- 
ferent from the estimate for the site diluted 3D Ising 
model (v = 0.684(5)). The same procedure was ap- 
plied to the magnetization evaluated at the temperature 
where the susceptibility is maximal. Judging the values 
of x 2 /d.o.f would lead to the result given in bold face in 
Table [j], but the effective exponent for the magnetization 
is clearly not yet stable. We therefore also considered 
the FSS behavior of higher (thermal) moments of the 
magnetization, (/j, n ) which should scale with a dimension 
n(3/v. The results for the first moments exhibit, how- 
ever, again much stronger corrections to scaling than we 
observed for \ or &k him, leading to our final estimate 
of 0/u = 0.65 ±0.05. . 

From the log-log plots of the three quantities in Fig. || 
one can clearly observe a crossover from a percolation- 
type behavior at small sizes, characterized by the expo- 
nents @ ~ijv ~ 2.05, l/v ~ 1.124, and j3/v ~ 0.45, 
towards a new regime at large sizes, which presumably 
corresponds to the random fixed point, with exponents 
as given above. The numerical evidence for this interpre- 
tation is quite striking, but with the present system sizes 
we can of course not completely rule out the possibility of 
corrections to scaling, in particular since for the 3D dis- 
ordered Ising model it is well known that such corrections 
at the random fixed point are strong (with a correction- 
to-scaling exponent around u> — 0.4). In order to inves- 
tigate this question for the 3D 4-state Potts model, we 
tried to fit the physical quantities to the standard ex- 
pression, e.g., a x L 1 l v \Y + t> x L~ u + •■•], including a sub- 
dominant correction-to-scaling term. Since 4-parameter 
non-linear fits are notoriously unstable, we performed lin- 
ear fits where the exponents are kept fixed and only the 
amplitudes are free parameters. In Fig. ^J, we show a 3D 
plot of the total x 2 for the susceptibility fits as a func- 
tion of j/u and u. We observe a clear, stretched valley 



3 



which confirms that 7/1/ is close to 1.5, but obviously this 
does not allow any reliable estimation of the correction- 
to-scaling exponent uj. The same procedure for 1/v gives 
qualitatively the same picture and confirms our previous 
estimate of 1/v = 1.33. For /3/v, on the other hand, 
the x 2 -landscape turns out to be very flat and extremely 
sensitive to the fit range. 




FIG. 4. Plot of the % 2 deduced from linear fits of 
Xma X (i) = a x L 7/V (l + b x L~ u ) in the range 25 < L < 96. 
The exponents are fixed parameters and the amplitudes are 
free. The base plane gives the ranges of variation of the expo- 
nents: 1.25 < 7/1/ < 1.75 and < w < 5. A cutoff at \ 2 = 50 
has been introduced for clarity of the figure. 

To conclude, from large-scale MC simulations of the 
3D 4-state bond diluted Potts model we obtained a) the 
phase diagram in the p — T plane in very good agreement 
with (rescaled) mean-field theory, b) the approximate lo- 
cation of the tricritical point around ptcp = 0.76(8), and 
c) at the dilution p = 0.56 clear evidence for softening of 
the rather strong first-order phase transition in the pure 
case towards a continuous transition with estimates for 
the critical exponents of v = 0.752(14), j/u = 1.50(2), 
7 = 1.13(4), fi/v = 0.65(5), and /3 = 0.49(5). These are 
clearly different from the values for both the disordered 
3D Ising [v = 0.684(5), 7/1/ = 1.963(5)] and the 3-state 
Potts [v = 0.690(5), 7/1/ = 1.922(4)] models. 
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